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Abstract 

We present an efficient approacli to study tlie carrier transport in graphene nanoribbon (GNR) devices using tlie 
non-equilibrium Green's function approacli (NEGF) based on tlie Dirac equation calibrated to the tight-binding n- 
bond model for graphene. The approach has the advantage of the computational efficiency of the Dirac equation 
and still captures sufficient quantitative details of the bandstructure from the tight-binding /r-bond model for 
graphene. We demonstrate how the exact self-energies due to the leads can be calculated in the NEGF-Dirac 
model. We apply our approach to GNR systems of different widths subjecting to different potential profiles to 
characterize their device physics. Specifically, the validity and accuracy of our approach will be demonstrated by 
benchmarking the density of states and transmissions characteristics with that of the more expensive transport 
calculations for the tight-binding /r-bond model. 
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1 Introduction 

Recent progress of graphene nanoribbon (GNR) fabrica- 
tion has demonstrated the possibihty of obtaining nano- 
scale width GNRs, which have been considered as one of 
the most promising active materials for next generation 
electronic devices due to their unique properties such as 
bandgap tunability via controlling of the GNR width or 
subjecting GNR to external electric/magnetic fields [1-5]. 
Device simulations play an important role in providing 
theoretical predictions of device physics and characteris- 
tics, as well as in the investigation of device performance, 
in order to guide the development of future device designs. 
Due to the nano-scale structures of GNRs, however, semi- 
classical treatments of carrier transport [6], which are the 
mainstay of microelectronics, are no longer valid. As a 
result, quantum transport formalism based on models 
incorporating detailed atomic structures, such as the ab- 
initio types [7-9], is needed for the proper simulation of 
these materials. Unfortunately, a full-fledge ab-initio ato- 
mistic model for carrier transport simulation is still very 
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computationally expensive and impractical even with the 
latest state-of-the-art computing resources. In this study, 
we therefore develop an efficient model in which a tight- 
binding Dirac equation (TBDE), calibrated with para- 
meters from the tight-binding 7r-bond model (TB-tt) 
[10-13], is used together with the non-equilibrium Green's 
function approach (NEGF) [14] to investigate transport 
properties of GNRs. We compare the density of states, 
DOS{E), and the transmission, T(E), of selected GNR 
devices for our TBDE model with that of the more expen- 
sive TB-TT model. Good agreement is found within the 
relevant energy range for flat band, Laplace and single 
barrier bias condition. We believe that our model and cali- 
brated data for a side selection of GNR widths presented 
in this article provided researchers in the quantum trans- 
port an accurate and practical framework to study the 
properties, particularly quantum transport in arbitrary bias 
conditions, of GNR-based devices. 

2 Model 

The Hamiltonian based on the Dirac equation [15,16] 
for graphene is given as: 
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H = 



where 



U{x) Vp{px — ipy) 
Vp{px + ipy) U[x) 



(1) 



-ihdu is the momentum for the direction 



^ = {x,y}, Vp is the Fermi velocity of graphene at the 
Dirac points (fixed at 10^ ms '^) and U{x) is the on-site 
potential. Due to the ID property of GNRs, the finite 
difference approach can be used along the transport 
direction (x) of GNRs and the Hamiltonian at each 
site «, and its backward and forward (h+) couplings 
with its neighbors (separated by a uniform spacing Iq) 
for a particular subband mode ky, can be written as: 



hn 



h- = [KV 



ihvpky 



Un 

ihvpky Un 
ihvp 
'2h 



0 1 
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(2) 



where Iq is the effective ID cell size as a result of the 
discretized Hamiltonian in (2). Figure la shows the 
schematics for real-space graphene and Figure lb the 
ID GNR model associated with (2). For an infinitely 
long GNR with uniform Uo, the Bloch waves solutions 
are valid and the dispersion relation, E{k„ky), for (2) is 



E[k„ky) = U0± 



hvp 



[kylo) +sin [kxh). 



(3) 



where for a fixed ky the positive and negative signs 
denoting the conduction and valence bands, respectively. 
In the absence of external potential {Uq = 0) and in the 

limit of large GNR width at which |fe| is small, (3) gives 

the linear dispersion for graphene E[k) = ±hvp |fe| . The 

energy bandgap of a certain width, and hence ky, is 
given by Eg = Ihvpky at k^ = 0. 

For non-equilibrium situations, we have to calculate 
the device retarded Green's function G(E) for a 




Cx 



►O 



Figure 1 Schematic representation of mapping of (a) a real- 
space two-dimensional GNR to (b) the one-dimensional Dirac 
Equation model with two degrees of freedom per effective cell 
of length £o 



particular energy E for the Hamiltonian in (2). Assum- 
ing the potential energies at the equilibrium source and 
drain are Us and Ud, respectively, and there are N lattice 
points in the device region, the G(E), of matrix size 2N 
X 2N, is given by G(£) = [EI2N - H - - i^d] ^ where 
the 'self-energies' and are associated with the 
effects of the semi-infinitely long source and drain [14]. 
Consider the self-energy of the drain (specified by the 
Hamiltonian of size 2M x 2M, where M is an arbi- 
trary number of lattice points with spacing Iq spanning 
the drain), defined in the NEGF framework [14] by 
Ed = t+C/(£)t_, where the drain Green's function, 
g{E) = {E- Hd)-^ , is also of the size 2M X 2M, and r_ 
= (r^^)^ is the coupling matrix (of size 2M x 2N) between 
the device and drain, which ends and starts at lattice 
points « = - 1 and 0, respectively. However, the only 
non-zero component of t+ is that of h± across the n = - 
1 and 0 interface, and hence only the 2x2 drain sur- 
face Green's function t/o,o, makes non-trivial contribu- 
tion to Sd, i.e., Od = h+Qofih- is the only non-zero 2x2 
submatrix, associated with lattice point « = - 1, of Z"^ (of 
size 2N X 2A0. Using the identity [Ehu - Hd)Q = hu 
for the drain region (« > 0), the system of equations for 
the dimensionless Green's function G can be written as 



Qo,o - h+Q\fl = h 



0 



5,1-1, 



+ OJ^°^Qn,o - KG, 



+yn+i,o 



0, n>l 



(4) 



(5) 



where ft)*°' = EI2 - ho is independent of sites inside the 
drain with uniform Ud- One can iteratively substitute 
Gn>o,o (second term) in (5) with the same in (4) so that 
after £ > 1 number of iterations, (4) and (5) can be 
rewritten as [17] 



(J^^o^Go.o =h+ a'-'^'Gi'.o 



(6) 



Qlhn.O = a^^' [^/2'(m-l),0 + ^/2'(m+l),o] ,m>l (7) 



where 



= K 



«w = /jw, £ > 1 



(8) 



(9) 



(10) 



«m=«(^-i)_2am = J2m(;^)«C0) (12) 



Chin et al. Nanoscale Research Letters 2012, 7:1 14 
http://www.nanoscalereslett.eom/content/7/l /II 4 



Page 3 of 7 



(13) 



The prefactor = g*«'o is such that is related to £ 
via (3). The integer m > 0 labels the surviving lattice 
points with spacing 2 '/q- The effects of the eliminated 
nodes after € number of iterations are taken into 
account in terms of "renormalized" couplings a*^^ and P 
(which happens to be equal in this model) and site 
energies (cy*^' at site index 2^m with m > 1 and co^^^ at 
m = 0, respectively). The symmetries of and h+ in (2) 
resulted in a^^\ co^^^ and co^^ each directly proportional 
to the "bare" energy for all i > 1, with their respec- 
tive coefficients A^^\ and ^^^^ as scalar functions 
dependent on X only. We show by induction that for all 
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Figure 2 The total computing time for calculating a series of 
Go.oiE) for all relevant modes in - 1 < £ < 1 eV with 0.001 eV 
spacing using analytic (o) and iterative (A) methods in the 
TBDE model for different GNR width. The iterative method takes 
about 40 X longer than that of analytic method. Included for 
comparison is the total time to calculate the corresponding surface 
Green's functions calculated using iterative method in the JB-n 
model (0). 



1-12 



1 + 



Z^-' (-1)^2, . 



(16) 



uniquely satisfy (11), (12), and (13). Since we are inter- 
ested in the retarded Green's function (i.e., £—>■£+ irj) 
for an infinitesimally small energy T] > 0, the condition 
imposed on the propagating waves is such that |A| = 1 - 
{lo/hvgjri < 1, where Vg = h'^ {dE/dkJ > 0 is the relevant 
group velocity [18,19]. Expanding in terms of A and tak- 
ing the Umit £ ^ oo, (14), (15), and (16) give A^°°^ = 0, 



n 



(-) 



(1 + A2)/(l - A2), and = 1/(1 - l2), respec- 
tively. The exact value of C/o,o> in the limit of € — > oo in 
(6), is now given by 



Go,o 



412 



12 - 1 



'0 

hvp 



2 r 



E-Ud 
ihvpky 



—ihvpky 
E-Ud 



(17) 



Similar argument can be applied at the source-channel 
interface where the analog source-side counterpart of 
Qo,o takes the same form as (17) with Us replacing Ud- 
Therefore, the only non-zero 2x2 submatrices for Z^^, 
d] are 



12 



ky 



ihVpky, 



-ihvpky E- U[s,d] 



(18) 



In the past, (6)-(13) are evaluated iteratively to calcu- 
late Go,o> and hence d] [13,17]. In this study, we 



have shown that (6)-(13) can be solved analytically for 
the Dirac form in (2) and that significant computational 
saving and accuracy can therefore be achieved by 
directly using (18) instead of numerically iterating (6)- 
(13). Figure 2 shows that the total computing time to 
calculate all the relevant modes of Go,o{E) for E e [-1,1] 
eV with spacing of 0.001 eV via analytical, i.e., (17), and 
iterative means, i.e., (6)-(13) for a range of GNR width 
on a typical duo core PC using MATLAB. The time 
needed to calculate Qo,o using the iterative method is 
about 40x larger than that of the analytic method over 
the entire range of the GNR width considered. In gen- 
eral, it is observed that the computing time increases 
with the GNR width for both analytical and iterative 
methods because the number of modes also increases 
with the width. (See Table 1.) Figure 2 also shows, as a 
comparison, the corresponding total computing time for 
calculating the all relevant surface Green's functions (via 
iterative method) for the same set of GNR width in TB- 
n model. This time is much larger than that of the 
TBDE, between about lOOx (at 1.1 nm width) and 455x 
(for 3.8 nm width) that of the analytic method of TBDE. 
Therefore the computational saving from using our ana- 
lytic results for the surface Green's function, (17), is 
compelling. The computing saving will be even more 
apparent in more realistic quantum transport calcula- 
tions in which the NEGF and Poisson equation are 
solved iteratively to achieve self-consistent solutions. 



Chin et al. Nanoscale Research Letters 2012, 7:114 
http://www.nanoscalereslett.eom/content/7/1/1 14 



Page 4 of 7 



Table 1 Results of best-fit Ig (for their respective subbands) to be used for our TBDE model for GNRs of different 
widths 



Family 
3p 


(eV) 


Subbands 
(eV) Ho (nm)] 


Family 
3p +1 


(eV) 


Subbands 
(eV) Ho (nm)] 


W12 


122 


0.612 [2.300], 
0.859 [1.860] 


WIO 


0.874 


0.437 [1.960], 
1.273 [2.712], 
1 .808 [1 .650] 


WIS 




0.477 [2.258], 
0.682 [1.917], 
1.654 [1.675], 
1.760 [3.150] 


W14 


0.675 


0.337 [2.002], 
0.966 [2.528], 
1 .446 [1 .725] 


W23 


U.oo 


0.331 [2.230], 
0.482 [1.974], 

1.208 [2.741], 

1 .209 [1 .800], 
1.831 [1.710] 


W18 


0.549 


0.275 [2.031], 
0.778 [2.442], 
1.201 [1.832], 
1.914 [3.150], 
1.963 [1.630] 



With G(£) now specified, the DOS(E), T(E), line 
charge density (pio) and total current (It) can be 
obtained, respectively [20], via 



DOS{E] 



In 



T(£) =Tr[r,GrdG'], 



(19) 



(20) 



/dEj 
^Diag[G(r/, + rrf/d)G-f], (21) 

SB " 



where T[,4] = i (^^[s.d] - ' f^.di^^ is the Fermi 

function at either the source or drain, Zsb denotes sum 
over the subbands, Diag[...] and Tr[...] denote the diago- 
nal and the trace of a square matrix, respectively. 

3 Results and discussions 

To incorporate the material details of GNR into the TB- 
n model, we first fit (3) of different GNR widths with 
that of the TB-tt model, which is widely used to calcu- 
late the bandstructures of GNR, for a flat potential (i.e., 
U = 0). Both real and imaginary parts of (3) are fitted 
for multiple subbands with different values of Iq for a 
particular GNR system. Figure 3 shows the comparisons 
of E{k) for the GNRs with width 1.0 nm and 1.4 nm, 
labeled as WIO and W14, respectively. At larger k, the E 
(k) calculated using (3) deviated from the that of the 
TB-TT model. This is expected as the TBDE model for 
GNR is most accurate near the Dirac points at small k 
[15]. Since we are interested in semiconductor proper- 
ties of GNRs, only the wide bandgap armchair GNRs 
(families with indices of w = 3p and 3p+l) [8,21] are 



considered here. The GNRs associated with w = 3/? + 2 
have Eg that are too small and are not considered here. 
Table 1 shows the best-fit Iq at different subbands for 
the m = 3p and m = 3p + I GNRs obtained under this 
study. With these calibrations, the adequate bandstruc- 
ture details based on TB-tt model can be incorporated 
in the TBDE model. Figure 4 compare the DOS(E) and 
T{E) for the same W12 and W14 systems using TBDE 
model (with the fitted-lo values from Table 1) and that 
of the TB-TT model. The very good agreements of results 
between the two models is a good first step to demon- 
strate the validity of the TBDE model in tackling quan- 
tum transport problems at which accurate T{E) and 
DOS{E) are the keys. 

To apply the NEGF-TBDE to more realistic transport 
situations, one needs to solve the NEGF-TBDE under 
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Figure 3 The £(/c) calculated from the TBDE matching that of 
the TB-7T with different best fit /q for different subbands for 

the (a) W12 and (b) W14 devices. Only conduction bands for f > 
0 are shown. The valance bands are symmetric about £ = 0. 
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Figure 4 The DOS(£) and T{E) of (a) W12 and (b) W14 devices 
with the best-fit /q (from E(k) calculations) at (7 = 0 agreeing to 
that of the TB-tt model. Both the DOS(Q and T{Ej are symmetric 
about f = 0. 



bias potentials. For a Laplace potential (with a bias of 
0.3 V), as shown in Figure 5a, the DOS{E) and T{E) for 
the W14 GNR are shown in Figure 5b, c, respectively. 
The corresponding TB-tt results and that of TBDE 
model with U = 0 are also included for reference. As 
shown in Figure 5a, the 0.3 V bias is achieved by shift- 
ing the conduction and valence bands upwards relative 
to those at the drain. As the highest valence band-edge 
(£v) (tit source) shifted up by 0.3 eV, the onset of DOS 
(E) for £ < 0 also creeped up into the original forbidden 
zone (with Zi = 0) by about 0.3 eV as indicated by arrow 
in Figure 5b. The positions of the DOS{E) associated 
with the higher subbands have also moved up the 
energy scale relative to those for U = 0. However, the 
onset of DOS(E) for £ > 0 has not been affected signifi- 
cantly by the Laplace setup because the lowest conduc- 
tion band-edge, which is at the drain, is still intact at 
E = EJ2. Although the forbidden zone for DOS{E) has 



channel 

Laplace 



W14 



E = 0 




E[eV] 



Figure 5 The 00S(£) and T{E) of a simple Laplace Potential (a) 

Schematic of a simple Laplace potential profile with a bias of 0.3 V 
across the GNR channel, (b) The resulting DOS(E) versus f with red 
arrow indicating the new addition of DOS(E) due to the upward 
movement of valence band-edge by 0.3 eV. (c) T(E) versus E. Results 
for U = 0 and that calculated by TB-tt are also included for 
comparison. 



narrowed as indicated in Figure 5b, the forbidden zone 
for T(E) has actually widen, as shown in Figure 5c, with 
the onset of non-zero T(E) for £ > 0 receding upwards 
by about 0.3 eV as indicated by the arrow, but 
unchanged T(E) for £ < 0. This is because from carriers 
are only unhindered source-to-drain only at £ >Eg/2 + 
0.3 eV and £ <Eg/2. The newly addition of DOS{E) at 
the source-side valence has no state of comparable £ to 
connect to in the channel and drain and hence does not 
contribute to £(£). 

Next, we subjected the W14 GNR to a rectangular bar- 
rier of 0.1 eV in the channel as shown in Figure 6a. The 
resulting DOS(E) and T(E) are shown in Figure 6b, c, 
respectively, with that of TB-tt model and U = 0 included 
for comparison. As expected, the onset of both DOS(E) at 
the conduction and valence ranges have not changed 
because the lowest E,. and highest £^, at -Eg/2, and Eg/2, 
respectively, have not been changed by the introduction of 
the barrier potential compared to that oi U = 0. However, 
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W14 



0.1 eV 



Eg = 0.675 eV E = 0 

' 0.1 eV 




-0.4-0.2 0.0 0.2 0.4 0.6 0.8 1.0 
E[eV] 

Figure 6 The DOS[E) and T{E) of a rectangular barrier, (a) 

Schematic of a rectangular barrier of height 0.1 eV across the W14 
GNR channel, (b) The DOS(E) with red arrow indicating region of 
reduced DOS(E) due to the introduction of the barrier at channel. 
The region near E = Eg/2 (as indicated) is magnified as inset with 
D05(E) in log-scale. Two discrete bound states, created by the 
inverted well at valence band-edge, as shown in (a), (c) The T(E) 
with red arrow indicates the receding T(E) away from Eg/2 due to 
the 0.1 eV barrier Results for (J = 0 and that calculated from TB-n 
are also included for comparison. 



it is observed that the magnitude of DOS(E) just above E = 
Eg/2 was reduced significantly due to the lost of states in 
the channel region dominated by the barrier. The inverted 
well of depth 0.1 eV at the channel valence band-edge is 
expected to accommodate some discrete bound states. 
However, the DOS(E) associated with them may be too 
sharp to be captured, or partially captured by the E grids 
being used. This expectation is borne out by the inset of 
Figure 6b, which shows the log-scale of the DOS(E) in the 
vicinity of £ = -Eg/2. Two discrete bound states, with the 
heights of their DOS(E) partially captured, are discernible 
within the inverted well energy range of within 0.1 eV 
above -Eg/2. As for T(E), the carriers are unhindered 
source-to-drain only for E >Eg/2 + 0.1 eV and E < -Eg/2 



eV and hence those boundaries marked the onset of T{E), 
as shown in Figure 6c. The bound states created by the 
inverted well in the channel region do not contribute to T 
(E) as there are no states of comparable energies both at 
the source and drain to connect to them. 

In both the Laplace and rectangular barrier potential 
profiles, the DOS{E) and T(E) for our TBDE model are 
in satisfactory agreement with that calculated from TB- 
n model within about 1.5 eV around the mid-gap. At 
higher energies, significant deviations in the DOS{E) and 
T(E) are consistent with the discrepancies we observed 
in £(/c) (as shown in Figure 3b), as discussed earlier. 
Nonetheless, these deviations are limited to the high- 
energy range that is of little relevance to the electron 
transport in GNR devices. Therefore, our TBDE 
approach is expected to be valid and as a practical and 
efficient alternative to TB-tt for studying carrier trans- 
port involving arbitrary self-consistent electrostatic 
potentials for device simulations [22,23]. 

4 Conclusion 

We developed a tight-binding Dirac equation for practi- 
cal and accurate numerical investigation of the electron 
transport in GNR devices. Based on our knowledge, this 
is the first time that the surface Green's function arises 
from applying the Dirac equation in NEGF framework is 
calculated exactly and hence can be used to achieve sig- 
nificant savings in NEGF calculations. The TBDE model 
is calibrated, with the appropriate parameters {vp = 10^ 
ms'^ and lo), to match the relevant bandstructure details 
as that of the TB-tt model, especially near the Dirac 
points. The best-fitted /q for a selected set of GNR 
widths are also presented for use. We show that the 
DOS{E) and T(E) calculated by our calibrated TBDE 
model can produce very good agreement with those that 
are calculated by the more expensive TB-tt model for 
the flat, Laplace, and rectangular barrier potentials. 
These cases validate the accuracy of the TBDE model 
and provided good confidence that the model can be 
used as a practical and accurate starting point for quan- 
tum transport of GNR-based devices where non-equili- 
brium and arbitrary electrostatic potentials are involved. 
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